Descriptive transcriptome analysis of tendon derived fibroblasts following in-vitro exposure to advanced glycation end products

Background Tendon pathologies affect a large portion of people with diabetes. This high rate of tendon pain, injury, and disease appears to manifest independent of well-controlled HbA1c and fasting blood glucose. Advanced glycation end products (AGEs) are elevated in the serum of those with diabetes. In vitro, AGEs severely impact tendon fibroblast proliferation and mitochondrial function. However, the extent that AGEs impact the tendon cell transcriptome has not been evaluated. Objective The purpose of this study was to investigate transcriptome-wide changes that occur to tendon-derived fibroblasts following treatment with AGEs. We propose to complete a descriptive approach to pathway profiling to broaden our mechanistic understanding of cell signaling events that may contribute to the development of tendon pathology. Methods Rat Achilles tendon fibroblasts were treated with glycolaldehyde-derived AGEs (200μg/ml) for 48 hours in normal glucose (5.5mM) conditions. In addition, total RNA was isolated, and the PolyA+ library was sequenced. Results We demonstrate that tendon fibroblasts treated with 200μg/ml of AGEs differentially express 2,159 gene targets compared to fibroblasts treated with an equal amount of BSA-Control. Additionally, we report in a descriptive and ranked fashion 21 implicated cell-signaling pathways. Conclusion Our findings suggest that AGEs disrupt the tendon fibroblast transcriptome on a large scale and that these pathways may contribute to the development and progression of diabetic tendinopathy. Specifically, pathways related to cell cycle progression and extracellular matrix remodeling were affected in our data set and may play a contributing role in the development of diabetic tendon complications.


Introduction
Tendon degeneration and impaired biomechanical function result in significant reductions in mobility and quality of life for the majority of the~30 million Americans living with diabetes, resulting in a substantial economic burden to individuals and society. Compounding the problem, human [1][2][3] and rodent [4] studies indicate that improving blood glucose levels does not normalize tendon properties in individuals with diabetes. Any new approach to enhance tendon health in people with diabetes is hindered by a poor understanding of the underlying etiology of tendon degeneration and impaired biomechanical properties [1,2,[5][6][7].
Our previous cell culture work implicated advanced glycation end-products (AGEs) as a potential mechanism driving tendon degeneration [8]. AGEs can form non-enzymatic crosslinks with collagen [9], a mechanism that has traditionally been the focus of tendon complications in persons with diabetes [10]. Yet, recent studies of tendons from humans with diabetes have found no evidence of greater collagen crosslinking than those without diabetes [1,11] and no relationship between tendon AGE content and tensile mechanics [11]. A less explored mechanism of AGE-mediated effects is the interaction of serum AGEs with AGE receptors (RAGE). AGEs accumulate in the serum of patients with diabetes [12][13][14] and our cell culture data suggest that AGEs can impact tendon cells. Specifically, treatment of cells with AGEs dose-dependently reduced cell proliferation and mitochondrial ATP production.
A thorough understanding of the cell signaling events contributing to the development of AGE-mediated diabetic tendinopathies will assist in exploring alternative areas of thought and developing therapeutic options to target this large patient population. Therefore, to better understand the effect of AGEs on tendon cells, we sought to characterize the alterations to the tendon fibroblast transcriptome following exposure to AGEs. Although many of these pathways have already been implicated with AGEs from analysis of non-tendon tissues, the primary goal of this study was to establish a descriptive and ranked evaluation of pathway disruptions that occur to tendon fibroblasts following an AGE insult.

Animal protocol
Animals utilized in this study were from a previous investigation [8]. The study was approved by the Purdue University Institutional Animal Care and Use Committee. All animals were cared for per the recommendations in the Guide for the Care and Use of Laboratory Animals. Five eight-week-old female Sprague-Dawley rats were purchased from Charles River Laboratories (Wilmington, MA) and maintained for an additional eight weeks. Rats were housed on a 12-hour light-dark cycle and provided standard rat chow and water ad libitum. At sixteen weeks (Final Weights: 256.43±5.19 g), rats were euthanized by decapitation after CO 2 inhalation. 0.2% type I collagenase, and incubated in a 37˚C shaking water bath for four hours. After digestion, the cell suspension was filtered through a 100μm mesh filter, pelleted by centrifugation, and resuspended in 5.5mM glucose DMEM containing 10% FBS, 1% sodium pyruvate (Sigma, St. Louis, MO), and 1% penicillin/streptomycin (Thermo Scientific, Waltham, MA). Samples were then plated in 100mm collagen-coated dishes. After reaching confluency, tendon fibroblasts were split and seeded (100,000 cells) in 100mm collagen-coated culture plates. Each donor animal's (n = 5) tendon fibroblasts were treated separately with 200μg/ml of BSA-Control or AGE-BSA for 48 hours for downstream paired DESeq2 analysis. Tendon fibroblasts treated at passages 2-4 were used for RNA isolation and RNA-sequencing (RNAseq).

Age preparation
Details on the preparation of AGEs have been reported previously [8,15]. Briefly, sterile filtered 30% BSA solution (Sigma, St. Louis, MO) was incubated with 70mM glycolaldehyde dimer (Sigma) in sterile PBS for three days at 37˚C. After incubation, the AGE product was dialyzed against sterile PBS for 24 hours at 4˚C using gamma-irradiated 10kDa cut-off cassettes (Thermo Scientific, Waltham, MA) to remove unreacted glycolaldehyde. Unmodified control BSA was prepared similarly, without the addition of glycolaldehyde dimer. Protein concentration was determined by BCA assay (Thermo Scientific) and absence of endotoxin (<0.25Eu/ml) was confirmed via the LAL gel-clot assay (GenScript, Piscataway, NJ). The extent of BSA modification was confirmed by fluorescence, absorbance, and loss of primary amines [15][16][17][18]. AGE-BSA and Control-BSA were diluted to 1mg/ml in PBS and fluorescent spectra and absorbance were recorded at 335nm excitation/420nm emission and 340nm, respectively (Molecular Devices, San Jose, CA). For determination of loss of primary amines AGE-BSA and Control-BSA were diluted to 0.2mg/ml in PBS. An equal volume of orthophthalaldehyde solution (Sigma) was added and fluorescent spectrum was recorded at 340nm excitation/455nm emission (Molecular Devices). Absorbance readings were completed to determine the extent of glycation. AGE-BSA showed increased glycation with absorbance readings of 0.682 AU compared to 0.01 AU for control BSA. AGE-BSA primary amine terminals underwent complete modification (-0.03% accessible amine terminals remaining), while control BSA retained 81.48% of accessible amine terminals. Negative values were interpreted as zero, and extent of modification was similar to previous reports [15].

RNA sequencing
Total RNA was isolated as previously described [8]. Briefly, RNA was isolated after BSA-Control or AGE-BSA treatment using the Direct-zol RNA Miniprep kit (Zymo Research, Irvine, CA). On-column DNase digestion was completed on all samples before elution of RNA. Total RNA from BSA-Control (n = 5) and AGE-BSA (n = 5) treated tendon fibroblasts was submitted to the Purdue University Genomics Core Facility (West Lafayette, IN) for PolyA + library construction. The integrity of input total RNA was assessed using a Bioanalyzer RNA Nano chip (Agilent 2100, Santa Clara, CA). Libraries from 500ng of input total RNA were constructed as directed by the Nugen Universal Plus mRNA-Seq + UDI kit (PN#9144-96), but the RNA fragmentation time was decreased from 8 minutes to 4 minutes. Final library products were subjected to a 0.7 Ampure:1 Sample ratio purification to reduce lower molecular weight

Bioinformatics
RNAseq raw data set quality and analysis was completed using Basepair software (New York, NY) pipelines. Reads were first aligned to the transcriptome derived from rn6 genome assembly using STAR with default parameters [19]. Next, read counts for each transcript were measured using featureCounts, and differentially expressed genes were determined using DESeq2 using a paired analysis [20,21]. An adjusted p-value cut-off of 0.05 (corrected for multiple hypotheses testing) was used. Finally, GSEA was performed on normalized gene expression counts, using gene permutations for calculating p-value. A log 2 fold change cut-off of 1.5 was enforced.

Descriptive pathway profiling
To preserve unbiased gene target selection and maintain a hypothesis-driven pathway selection, GeneGlobe (Qiagen, Hilden, Germany) pathway database was utilized to complete a descriptive approach to pathway analysis. We generated heat maps based on GeneGlobe RT 2 profiler arrays independent of whether those gene targets were significantly altered in our dataset. Gene targets in the RT 2 profilers but not in our dataset were excluded from heat maps. The percentage of significantly altered genes, both increased and decreased, was calculated based on the number of total genes included in each pathway's respective heat map to rank the most implicated pathways. This systematic approach was employed to maintain an objective view of the global data.

Pathway analysis
RNAseq data were imported into Ingenuity Pathway Analysis (IPA, Qiagen) to determine select pathways and biological functions that were altered in response to AGE-BSA treatment.

Results and discussion
Overview A total of 2,159 genes within our data set met the criteria of q<0.05 and fold change of greater than or less than 1.5 (log 2 fold change greater than or less than 0.584). One thousand forty-six genes were significantly increased, and 1,113 were significantly decreased (Fig 1).

Most affected gene targets
The top ten increased, and the top ten decreased gene targets within our data set were identified based on our log 2 fold change and adjusted p-value thresholds. The top ten increased gene targets in order of highest to lowest positive log 2 fold change were Cyp1a1, Pipox, Btc, Slc22a14, Tbxas1, Itgb2, Slc13a3, Cldn1, Ncf1, and Tnfrsf17 (Table 1). The top ten decreased gene targets in order of highest to lowest negative log 2 fold change were Pimreg, Pmch, E2f7, Pbk, Parpbp, Ube2c, Troap, Cenpf, Cldn23, and Ccnb2 (Table 1).

PLOS ONE
severe limitations to tendon fibroblast proliferative capacity and mitochondrial function while increasing mitochondrial DNA content [8]. We have followed up on these previous findings by completing a descriptive transcriptome profile of Achilles tendon-derived fibroblasts following AGE exposure. The goal of this study was to identify and rank pathways that were most implicated following AGE exposure, thus providing a more precise mechanistic exploration of AGE-mediated effects on tendon-derived cells.
Using a clinically-relevant concentration of AGEs [12,22], we have previously demonstrated incorporation of synthetic nucleoside 5-ethynyl-2´-deoxyuridine (EdU) in tendonderived fibroblasts to be~3% following AGE-BSA (200μg/ml) exposure as compared to~53% in the BSA-Control exposed group, which proliferate normally [8]. Further, we noted a  reduction in proliferative gene markers, Mybl2 and Pcna, and reduced absorbance values of cytostatic MTT with AGE-BSA treatment in tendon fibroblasts. Our RNAseq data corroborated our previous findings of reduced Mybl2 and Pcna gene expression and revealed several additional genes responsible for cell cycle progression to be significantly impacted (Fig 2). In fact, our transcriptome analysis revealed that genes associated with the cell cycle are the most impacted by AGE treatment (Fig 2 and Table 2). Tendon fibroblast proliferation is vital for tendon development and adaptation [23,24]. The inability of tenocytes to proliferate in the presence of AGEs could precipitate the development of tendon degeneration by limiting adaptations to loading [25]. Tendon healing requires a phase of increased cellular proliferation [23,24,26], thus AGE-induced limitations in cell proliferation could contribute to delayed in healing noted in those with diabetes [27][28][29]. In fact, would healing was identified as one of the top 10 GeneGlobe Pathways impacted by AGE treatment (Table 2 and Fig 10).
Gene targets associated with ECM maintenance and remodeling were also dramatically affected in our dataset (Fig 3). The ECM is vital to tendon tissue health and serves several vital functions, including cell adhesion, communication, and differentiation. Additionally, the ECM provides structural and biochemical support to the surrounding resident cell population. The tendon ECM consists primarily of type I and type III collagen fibers surrounded by proteoglycans that assist collagen fibrils' assembly and stability [30]. A precise and linear arrangement of collagen fibrils is vital to tissue integrity and, therefore, mechanical function [31]. The inclusion of multiple collagen isoforms allows the ECM to specialize and adapt to specific mechanical loading and functional responses [32]. For instance, type I collagen (Col1a1) is a stronger collagen isoform. In contrast, type III collagen (Col3a1) is weaker and generally upregulated in the early stages of tissue remodeling following exercise or during the initial stages of healing [33,34]. Col3a1 can provide temporary tensile strength to the tissue assembly until it is later replaced by stronger Col1a1 [35]. Although Col1a1 mRNA was unaffected in our RNAseq data set, Col3a1 mRNA expression was increased (Fig 3). Similarly, our previous report indicated Col3a1 mRNA expression increased with 50μg/ml and 100μg/ml AGE exposure compared to an equal dose of BSA-Control [8]. This increase in Col3a1 mRNA expression is likely in response to the AGE insult and an attempt to maintain the ECM environment.
Further, the most abundant tendon proteoglycan gene expression of decorin (Dcn) increased in our RNAseq data set (Fig 3). Dcn aids in the maintenance and regulation of collagen fibril structure and resident fibroblast proliferation [31]. As a critical regulator in matrix assembly, loss of Dcn would likely prove to be unfavorable to the strength of the tendon assembly, which would decrease the tissue's ability to withstand sudden strain [31]. Our observed increase in Dcn gene expression may be a compensatory response that results in response to the AGE insult. However, impacts to Dcn content and gene expression would need to be externally validated in a whole diabetic tendon.
Lysine and hydroxylysine are found within the collagen amino acid sequence and play an essential role in cross-link formation. Oxidation of lysine and hydroxylysine by lysyl oxidase (Lox) forms cross-links within collagen fibrils, contributing to tissue integrity by increasing tensile strength and stabilizing the collagen fibril assembly. Strength and stability of the tissue assembly are essential, especially given the high contractile forces tendons are responsible for transmitting from muscle to bone. Our dataset revealed Lox gene expression to be significantly reduced following AGE exposure (Fig 3). If reduced mRNA expression of Lox coincides with reduced enzymatic cross-link formation, AGEs may contribute to a weakened tendon assembly due to loss of enzymatic cross-links between adjacent collagen fibrils. Tendons of diabetic animals generally have a reduced load to failure capacity, which may be a result of greater tissue degeneration at the macroscopic level [4,11,28,36]. More work is needed to determine the impact of AGEs on the whole tendon fibril assembly. Remodeling of the ECM is primarily regulated by enzymes known as matrix metalloproteinases (MMPs), which are responsible for the degradation portion of ECM remodeling. Collagenases such as MMP-1 and MMP-13 cleave type I collagen molecules in the ECM. Similarly, gelatinases, such as MMP-2 and MMP-9, degrade collagen isoforms in the ECM. MMPs are transcribed and translated as proenzymes and then secreted into the ECM, where they are activated through proteolytic cleavage of the N-terminal. Although MMP activity is degenerative, it facilitates ECM remodeling and tendon tissue adaptation. In turn, MMP activity can be reversibly inhibited by a group of enzymes known as tissue inhibitors of metalloproteinases (TIMPs). TIMPs play an essential role in ECM remodeling by limiting MMP activity and preventing excessive degradation. Counter-regulation via TIMP activity tightly regulates the breakdown and synthesis of collagen in response to external stresses, such as mechanical loading. Loss of ECM regulation, such as favoring degradation over synthesis, could alter the ECM responses to damage the tissue assembly. It is no surprise that the dysregulation of degenerative enzymes, such as MMPs, has been thought to play an essential role in developing tendon pathology in diabetes as overexpression of MMPs may favor ECM degradation [37]. Similarly, if inhibitory TIMPs are less expressed, the environment may also favor degradation by allowing MMPs to act on the ECM for a more extended period. Previous reports have indicated that AGEs increase MMP -2, -3, -9, and -13 secretion and expression in chondrocytes with 100μg/ml of AGEs [38,39]. Further, mRNA expression of MMP -1, -3, and -13 in porcine chondrocytes was increased with 100μg/ml of AGE exposure [40]. Our previous work in Achilles tendon-derived fibroblasts demonstrated an increase in MMP -2 and -3 but no significant changes to MMP-9 and -13 [8]. Our RNAseq analysis confirmed MMP -2 and -3 to be elevated, along with MMP -15 and -17. However, we did not observe any changes to TIMP -1, -2, -3, or -4 in our RNAseq dataset, suggesting that MMPs may be exerting their function in an unorganized fashion that would favor a degenerative ECM environment. MMP gene expression data is limited in scope as it does not account for ECM secretion and N-terminal cleavage. However, the large impact that AGE exposure has on the dysregulation of ECM-related gene expression is further evidence that elevated serum AGEs may be contributing to the development of connective tissue pathology in diabetic populations (Fig 3). Delayed and abnormal healing is a common complication of types I and II diabetes [27,41]. Not only does it appear that diabetic patients are at risk of developing tendon tears, but healing post-repair is also impaired [42][43][44]. Interestingly, transforming growth factor (TGF) β1 expression was significantly reduced in our RNAseq data (Fig 7). In addition to TGFβ1 being one of the affected genes in the wound-healing pathway (Fig 10), the GeneGlobe TGFβ signaling pathway was also strongly influenced by AGE treatment (Table 2 and Fig 7). TGFβ is a critical factor in fibrosis and modulation of ECM homeostasis [45]. It has previously been demonstrated that TGFβ levels are significantly reduced in diseased human rotator cuff tendon samples [45]. In addition, TGFβ is known to modulate inflammatory responses by influencing fibroblast recruitment and stimulating collagen production [46,47]. Inconsistent with known effect of TGFβ on collagen production [46,47], Col1a1 was unchanged in our RNAseq dataset, and Col3a1 was increased (Fig 3). However, mRNA expression of Col5a1, Col5a2, and Col5a3 expression was significantly reduced in our RNAseq dataset. Type V is a fibrillar collagen isoform found less abundantly in a tendon but exists to provide support to tissues that do contain high levels of type V collagen isoforms [48]. While the wound healing GeneGlobe pathway was not as affected as other pathways, it is still likely that these gene targets contribute in some manner to the delayed healing response that is commonly observed following tendon injury in diabetic patients.

Conclusions
Several studies have shown that the risk of developing tendinopathy is greater in those with diabetes mellitus [42][43][44]49]. Our new data highlights cell-signaling pathways that may assist with expanding our understanding of diabetic tendon pathology and failed healing responses. While our discussion is limited in scope, and we provide only transcriptome data, the purpose of this study was to complete a descriptive profile of the AGE insult to tendon fibroblasts. This work is the first data set to utilize RNAseq methodology to study the tendon fibroblast transcriptome following AGE exposure. These data will be helpful for further elucidation of the diabetic tendon disease process.